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We study the shape of the tail of a heap of granular mate- 
rial. A simple theoretical argument shows that the tail adds 
a logarithmic correction to the slope given by the angle of re- 
pose. This expression is in good agreement with experiments. 
We present a cellular automaton that contains gravity, dissi- 
pation and surface roughness and its simulation also gives the 
predicted shape. 
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In this letter we discuss the shape of static heaps of 
granular material. In fact many typical phenomena ob- 
served in granular materials have been studied in the last 
years |l] ||. In order to understand these phenomena 
simple models as well as computer intensive simulation 
techniques jij have been very useful. 

As everybody knows real sand piles are almost perfect 
cones with a well defined angle of repose which depends 
on gravity and on the characteristics of the material. 
Watching carefully however, one notices the existence 
of curved tails at the bottom of the heaps to which in 
fact not much attention has been paid in the literature 
up to now. In the following we shall study the shape of 
these tails in more detail. 

We will present a simple experiment to obtain the pro- 
file of the tails in two dimensions. Then, we deduce an 
analytical expression for their shape which agrees with 
our experimental profiles. This expression has two pa- 
rameters that describe the material properties. Finally, 
we test the theory with data obtained from computer 
simulations using cellular automata (CA). 

Several CA models have been proposed in the past to 
study the behaviour of granular matter ||] . They in- 
clude dissipation which is the most important ingredient 
to capture the peculiarities of granular media. One such 
automaton is a lattice gas formulation by Peng et al. ||] 
with rest particles and inelastic collision rules. Using 
this model we have simulated two dimensional heaps in 
a quasi-static regime and compared the resulting profiles 
with those predicted by our theory. 

We built sand piles in an easy and inexpensive exper- 
iment. Grains were poured at a slow rate of about ten 
particles per second from the top into the center of a rect- 
angular cell made of two parallel vertical glass plates of 
size 30 x 25 cm 2 separated by a fixed distance of 2 mm. 
As granular materials we used lead spheres, sugar and 



polenta. The grain diameter was 2 mm in the first case 
and about 0.5 mm in the other two cases. In fig. 1(a) 
we see a digizited image of the experiment showing the 
essentially two-dimensional heaps that we obtain at dif- 
ferent times. We have studied the profiles of the heaps by 
recording digitized pictures taken by a VHS video cam- 
era at different stages of the evolution. The resolution of 
digitized images was 40 pixels/cm. 

The heaps are grown in steady state, in the sense that 
only a few grains move simultaneously being in the con- 
stant velocity regime || . In this regime particles do not 
accumulate at the top of the heap and the formation of 
big avalanches which can disturb the formation of smooth 
profiles is not present. Therefore profiles do not depend 
on the flow rate at the top. They only depend on the 
type of grains used and gravity. 

In Fig. 1(a) we see the typical situation of two flat sur- 
faces that define the angle of repose. We will, however, 
focus our attention on the small curved tails that can be 
observed at the base of the heaps and obtain the analyt- 
ical expression of their shape using a simple argument. 

Taking pictures at different times we discovered that 
it is possible to superpose the profiles of the tails by a 
simple horizontal shift (see Fig. 1(a)). It seems that the 
heap grows by putting layers of particles over one an- 
other. Each individual layer grows upwards from the 
bottom to the top by stopping particles which are mov- 
ing down along the surface of the heap. 

At the top of each incomplete layer there is a kink 
(kinks are marked by arrows in Fig. 1(b) ^ The presence 
of such kinks reduces the slope of the surface away from 
the angle of repose of the material. Let us describe the 
surface by h(x) where h is the height and x the corre- 
sponding horizontal displacement. We choose the origin 
at the center of the base of the heap, i.e. h(0) = h m 
where h m is the height of the heap at the top, and con- 
sider only x > since the heap is symmetric with re- 
spect to the origin. In the absence of kinks one would 
have h = h m — 72; where 7 = tan 9 and is the angle 
of repose. The presence of each kink increases this ideal 
value of x by a certain value I which typically is of the 



*In the case of irregulary shaped particles or particles of 
different size one has asperities of different size on the surface 
of the heap and the larger ones effectively act as the kinks 
discussed here. 
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size of a grain. If p is the number of kinks per unit length 
in the vertical direction we can express the slope of the 
surface as a function of p as: 
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Particles falling down along the surface collide with the 
kinks and can be accumulated on top of them with some 
rate r or continue to move downwards to the next layer. 
r depends on I and on the properties of the grains (shape, 
roughness, coefficient of restitution). Let us call the 
flux of particles pulled down by gravity along the surface. 
Since the experiments were performed in steady state one 
has the relation 4? = r$p. The fact that the heap grows 
by a shift of the profile in the horizontal direction means 
that the aggregation rate of particles is independent of 
h. Since the number of particles aggregated per vertical 
unit length is the variation of the vertical flux we have 
d&/dh = B, where B is a constant and since $(0) = Owe 
obtain $ = Bh. Putting everything together the slope of 
the surface as a function of p is given by: 
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After integration we obtain the final expression, 

x = h l e >- 9—r- 

7 h 



(2) 



(3) 



which shows a logarithmic correction to the angle of re- 
pose. 7 = tan 9 is the classical dynamic friction coef- 
ficient and l e = - describes the typical extra horizon- 
tal displacement a particle must undergo before it sticks. 
Both constants characterize the granular material. In 
Figure 1(c) we show the average over 10 different profiles 
representing the deviation from the straight profile. Er- 
ror bars are the standard deviations for each value of h. 
The surface fit is for 7 = 0.98 and l e = 1.5mm. It is very 
interesting to note that for the three materials used here 
we find consistence with r « 1/3. 

In the following we will check our formula using numer- 
ical simulations. We consider a Lattice Gas Automaton 
(LGA) at integer time steps t = 0, 1, 2, ... with particles 
located at the sites of a two dimensional triangular lattice 
of size L. Gravity acts downward in the vertical direction 
and forms an angle of 30° with the closest lattice axis. 
At each site there are seven bit variables which refer to 
the velocities Vi(i — 0,1, 2, ...6). Here Vi(i = 1, ...6) are 
the nearest neighboring (NN) lattice vectors and vq = 
refers to the rest state (zero velocity). Each state can be 
either occupied by a single particle or empty Therefore, 
the number of particles per site has a maximal value of 
seven and a minimal value of zero. 

The time evolution of the LGA consist of a collision 
step and a propagation step. In the collision step parti- 
cles can change their velocities due to collisions and in the 



propagation step particles move in the direction of their 
velocities to the NN sites where they collide again. The 
system is updated in parallel. Only the collisions speci- 
fied in Figure 2 can deviate the trajectories of particles 
from straight lines with probabilities depending on the 
dissipation parameter p. If p is not zero, energy can be 
dissipated during the collision. This is a crucial property 
of granular materials and yields among other an instabil- 
ity towards cluster formation Jio[ . 

The two collisions shown on the lower part of Fig. 2 
temporarily allow more than one particle on a site. How- 
ever, immediately after the collision step, the extra rest 
particles randomly hop to NN sites until they find a site 
with no rest particle and there they stop. 

We incorporate the driving force, namely gravity g, 
by setting a rest particle into motion with probability 
g /2 along one of the two lattice directions which form an 
angle of 30° with the direction of gravity, however, only 
if the site below is empty at that time. Rest particles 
that are already on the heap can only be accelerated by 
gravity if at least one of the two NN sites in the direction 
of gravity does not yet belong to the heap. 

During the evolution, we add particles at a fixed rate 
from one site at the top to the system. Gravity moves 
these particles down until they collide with the hard wall 
at the bottom. A particle colliding with the bottom is ei- 
ther reflected with probability 1 — p in the specular direc- 
tion or stopped with probability p losing its momentum. 
In the second case the resulting rest particle belongs to 
the growing heap, and we label that particle in order to 
store the information that it is sustained by the bottom 
plate. Every particle colliding with these particles looses 
its momentum and is aggregated to the heap after the 
redistribution process described above has taken place. 

We perform simulations for systems of size L = 512 
typically iterating 3 x 10 5 time steps. Most of the com- 
puter work was performed on a CM-5. Like in the exper- 
iment we add few particles (one every 8 time steps) to 
be in a quasi-static regime. In that case the profile does 
not depend on the history of the system. 

For all parameter values the profiles have an angle of 
repose of 60° corresponding to 7 = y/H which is deter- 
mined by the geometry of the underlying lattice. The tail 
shows the presence of kinks which reduce the slope. The 
different gray levels correspond to different time steps. 
One notes that like in the experiment the surfaces of the 
heap are just horizontally displaced in time. 

For this model the parameters in equation (3) are easily 
determined. On one hand we use as unit length the dis- 
tance between NN sites on the underlying lattice which 
we choose to be unity giving 1 = 1. On the other hand 
r = 1/3 because every particle colliding with a kink is 
aggregated to the heap and then redistributed randomly 
to one of the three empty NN sites. Only one of these 
sites particles will aggregate on top of a layer. 

In Figure 3 we have plotted an average of 16 indepen- 
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dent profiles in the same manner as in Fig. 1(b). The 
continous line is the shape resulting from equation (3) 
with 7 = \/3 and l e = 3. Our formula is in excellent 
quantitative agreement with the simulation. 

In order to include into the CA the sticking of parti- 
cles typical for wet sand we introduce a new parameter, 
namely the probability ij that one particle stopped and 
aggregated after a collision with the heap is "blocked". 
Gravity can only move a "blocked" particle if none of 
the two NN sites below belongs to the heap and then the 
particle is no longer blocked. Introducing these rules the 
slope of the surfaces of the heaps can be larger than 60° 
because each blocked particle can support rest particles. 

Averaged profiles for r\ = 0.002 and 0.004 are shown 
in Fig. 4. We observe that now there is no well defined 
angle of repose like in the case of dry sand. The profile 
can be calculated using very similar arguments to those 
used before. 

Let po be the vertical density of blocked particles on 
the surface. The relation between the slope of the profile 
and po is then 
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because the presence of each blocked particle on the sur- 
face reduces the value of x by I = 1. 

It is easy to find a relation between p and h. On one 
hand po oc ??<!>. On the other hand $ oc h due to the 
translational invariance which we have observed experi- 
mentally so that po = c'rjh. One therefore obtains, 



dh 

dx 
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1 + V3(^ - c' V h) 
which after integration gives 
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where c = ^dr\ is a constant. In Fig. 4 we also show 
the fits obtained from this expression using c' = 0.39 and 
r = 0.35. We have checked that po oc h in agreement with 
the argument for h < 100. If h becomes larger than 100 
the density of blocked particles saturates to po ~ 0.28. 

We have studied experimentally, theoretically and nu- 
merically with a cellular automaton the tail of sand piles. 
The theoretical argument giving the shape of the tail is 
based on the experimentally observed translational in- 
variance and uses mass conservation. One finds that 
the tail constitutes a logarithmic correction to the naive 
straight slope given by the angle of repose. The analytic 
expression fits very well to the experimental shape and 
to the two-dimensional heaps obtained with a cellular 
automaton. 



Acknowledgments 

We are indebted to S. Roux for fruitful discussions, to 
P.Petitjeans and D.Hoang for help on the experimental 
work and to H. Puhl for much help on the computer. JJA 
thanks the CNCPST for a generous grant of computer 
time on the CM-5 and is grateful for an Ayuda parcial 
(PB91-0709) and a postdoctoral grant from DGICYT. 



[1] H.M. Jaeger and S.R. Nagel, Science 255, 1523 (1992). 
[2] A. Hansen and D. Bideau (eds.), Disorder and Granular 

Media (North Holland, Amsterdam, 1992). 
[3] A. Mehta (ed.), Granular Matter (Springer, Heidelberg, 

1994). 

[4] H.J. Herrmann, III Granada School, Lecture Notes in 
Physics, vol. 448, Springer Verlag (1994). 

[5] G.W. Baxter and R.P. Behringer, Phys. Rev A 42, 1017 
(1990); Physica D 51, 465 (1991). 

[6] A. Karolyi, and J. Kertesz, Proceedings of the 6th Joint 
EPS- APS International Conference on Physics Comput- 
ing, (1994); S. Vollmar and H.J.Herrmann, preprint. 

[7] D. Desirable and J. Martinez, in Powder and Grains, ed. 
C. Thornton (Balkema, Rotterdam, 1993), p. 345; 

[8] G. Peng and H.J.Herrmann, Phys.Rev. E 48, R1796 
(1994) and preprint. 

[9] G.H. Ristow, F.-X. Riguidel and D. Bideau, J. Physique 
I, 4, 1161 (1994). 
[10] I. Goldhirsch and Z. Zanetti, Phys. Rev. Lett. 70, 1619 
(1993). 



3 



Figure captions 



[Fig. la, b, c] (a) Digitized image of a heap of the 
polenta. The diameter of the grains is 
about 0.5 mm. The height of the heap 
is 16 cm. Different gray levels show the 
pile at different stages of growth. The 
superposed continous lines in both fig- 
ures are fits obtained from eq. (3) by 
taking the values 7 = 0.98 and l e — 
1.5 mm. (b) Tail of a pile made of lead 
spheres with a diameter of 2 mm. One 
can observe parallel layers terminating 
in horizontal kinks (marked by arrows), 
(c) Deviation of the shape from the 
straight profile for the same material. 
Data are averages obtained from 10 dif- 
ferent samples. Error bars represent 
the standard deviations of the averaged 
values. Continuous line represent the 
same fit of Fig. (a) 

[Fig. 2] Collision rules of the cellular automa- 

ton. Arrows represent moving particles 
and full dots stand for rest particles. 
The number next to each configuration 
is the probability for that transition. 

[Fig. 3] Deviations from from the straight line 

for a heap of 80000 particles obtained 
in a quasi-static regime adding one par- 
ticle every 8 time steps. Data are av- 
erages obtained from 16 independent 
samples. The fit is the same as in 
fig 3(a) 

[Fig. 4] Shapes obtained from 16 independent 

samples for a roughness parameter r\ = 
0.002 (circles) and 77 = 0.005 (squares). 
Errors are of the same size of symbols. 
Continuous lines are fits obtained from 
eq. (6) with d = 0.39 and r = 0.35. 
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